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ABSTRACT 

A rece nt approach t o sim ulating localized feedback from active galactic nuclei 



(AGN) by Power et al. (20111 uses an accretion disc particle to represent both the 
black hole and its accretion disc. We have extrapolated and adapted this approach to 
simulations of Milky Way-sized galaxy mergers containing black holes and explored the 
impact of the various parameters in this model as well as its resolution dependence. The 
two key parameters in the model are an effective accretion radius, which determines 
the radius within which gas particles are added to the accretion disc, and a viscous 
time-scale which determines how long it takes for material in the accretion disc to 
accrete on to the black hole itself. We find that there is a limited range of permitted 
accretion radii and viscous time-scales, with unphysical results produced outside this 
range. For permitted model parameters, the nuclear regions of simulations with the 
same resolution follow similar evolutionary paths, producing final black hole masses 
that are consistent within a factor of two. When comparing the resolution dependence 
of the model, there is a trend towards higher resolution producing slightly lower mass 
black holes, but values for the two resolutions studied again agree within a factor of 
two. We also compare these results to two other AGN feedback algorithms found in 
the literature. While the evolution of the systems vary, most notably the intermediate 
total black hole mass, the final black hole masses differ by less than a factor of five 
amongst all of our models, and the remnants exhibit similar structural parameters. 
The implication of this accretion model is that, unlike most accretion algorithms, a 
decoupling of the accretion rate on to the black hole and the local gas properties is 
permitted and obtained; this allows for black hole growth even after feedback has 
prevented additional accretion events on to the disc. 

Key words: black hole physics - galaxies: interactions - galaxies: active - methods: 
numerical 



1 INTRODUCTION 

Observational evidence suggests that supermassive black 
holes exist at the centre of all galaxies with stellar spheroids 



(e.g. Kormendy & Richstone 1995 Ferrarese & Merritt 



2000 


). Further evidence, such as the Mbh-o" relationship 


(e.g. 


Silk & Rees|1998 


Ferrarese & Merritt|2000| |Gebhardt 


ct al. 20001 ITremaine et al.|2002||King'2003; Giiltckin et al. 


2009 


) and the MoH-Mbuige relationship (Magorrian et al. 


1998 


McLure & Dunlop|2002 Marconi & Hunt|2003l, sug- 



gests that the black hole and spheroid evolution are related. 
Moreover, galaxies are relatively smaller today than might 
be naively predicted from the hierarchical model of galaxy 
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formation, indicating that some mechanism has limited the 
growth of these galaxies. One favoured explanation is that 
during mergers gas from the merger fuels both star forma- 
tion and AGN activity (e.g. [Sanders et al.|19 88', Scannapieco| 
|et al.|[2005[ ). The latter is likely a self-regulatcd process: 
outflows from the black hole following a strong accretion 
event interact with the surrounding gas, inhibiting further 
accretion events, and hence limiting black hole growth (e.g. 
Silk & Rees"1998' |Fabian||1999l [Scannapieco fc Oh|[2n04l 
King,2005, ,King, 2010[ ). Ultimately, the feedback from the 
increased AGN activity blows away some or possibly all the 
gas, truncating the star formation and leading to an ellipti- 
cal galaxy (e.g. [Springel et al.|2005a 



Using various models, AGN feedback has been imple- 



mented in many numerical simulations (e.g. Springel et al 
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2005b| 


Thacker et al.||2006| |Okamoto et al.||2008| |Booth &| 


Schaye 


2009 


Kurosawa et al.|2009 


Debuhr et al.|2011 


1. The 



goals of this research have been varied, from reproducing the 
observed relationships between the black hole and spheroid, 
to other factors impacted by AGN activity, such as galaxy 



cluster properties (e.g. Puchwein et al. 2008 \ or even the 
impact on CMB foregrounds (e.g. Scannapieco et al.|2008 |. 

In all these simulations, the accretion rate of gas on to the 
black hole is dependent on extrapolating the macroscopic 
gas properties to the microscopic scale around the black 



hole. This is true for viscous accretion (e.g. Debuhr et al 
2011 Hopkins fc Quataert||20lT l, drag accretion (Okamoto 



et al. 2008), and the numerical form of Bondi accretion rate 
( |Bondi||1952| ), 



Mb 



(1) 



where a is a numerical parameter used to account for the 
limited dynamic range in the simulations, p and Cs are the 
gas density and sound speed around the black hole, Viai is 
the relative velocity between the gas and the black hole, and 
Mbh is the mass of the black hole. This accretion rate has 
been used in merger (e.g. Springel et al.||2005b Di Matteo 



et al. 2005 1 and in cosmological (e.g. Sijacki et al. 2007 
Booth fc Schaye|2009 \ simulations 



Most algorithms also limit the mass accretion to be no 
greater than the Eddington accretion rate. 



Me 



(2) 



where rup is the proton mass, ctt is the Thomson cross sec- 
tion, and Er is the radiative efficiency (i.e. the mass-to-energy 



conversion efficiency); we set — 0.1 (Shakura & Sunyaev 
1973). This rate differs from the previously considered ac- 



cretion rates in that it depends only on the black hole mass, 
and not the local gas properties. 

One of the assumptions frequently used in these models 
is that any accreted gas is immediately transferred to the 
black hole. This is clearly an unrealistic assumption as ma- 
terial must first shed its angular momentum. Indeed some 
models argue that understanding this process is one of the 
key factors in modelling AGN feedback; for example, the vis- 
cous accretion rate in Debuhr et al. (2011 1 accounts for an- 



gular momentum of the gas. However, the gas in this model 
is nonetheless instantly accreted on to the black hole. In re- 
ality the gas would be expected to settle on to a circular 
orbit of radius -Rcirc, which is set by the angular momentum 
of the gas and whatever processes cause it to shed angular 
momentum while collapsing (cf. Hobbs et al.|2011 1, resulting 
in an accretion disc around the black hole. The gas with the 
lowest angular momentum would then travel through the 
disc and eventually be accreted on to the black hole (e.g. 
King|2010[ ). 



Motivated by these ideas. Power et al. (2011 1 (hereafter 
PNKll) have implemented a two-stage accretion algorithm 
using an accretion disc particle (ADP), which models both 
the black hole and accretion disc. In the first stage, nearby 
gas, the bulk of which will have shed large amounts of an- 
gular momentum to reach this radius, is accreted on to the 
accretion disc; this accretion rate is dependent only on the 
relative positions of the black hole and the gas particle. In 
the second stage, the gas is accreted from the accretion disc 



on to the black hole. This method incorporates the delay 
between the time gas is accreted on to the disc and when 
it is finally accreted on to the black hole, and at the same 
time decouples the black hole's accretion rate from the in- 
stantaneous gas properties around the black hole. 

While the ADP model implemented in [Power et al.| 



(20111 relies upon resolving scales far smaller than those 
that can be resolved in merger simulations, the two stage 
approach, and perhaps more specifically the incorporation 
of a delay period before accreting on to the black hole, is an 
issue worthy of investigation in merger simulations which 
can have time-steps in the few thousand year range. But 
this consideration highlights the fact that resolution depen- 
dence must be considered. For example, at very low mass 
resolution and the associated low time resolution, concerns 
about accretion delays are likely less significant as the ratio 
of time-step and delay time get closer to unity. Therefore in- 
vestigating the precise resolution dependence is important. 
We also emphasize that the ADP model needs additional 
features to be implemented in a merger simulation, but we 
have drawn on prior work and use methods that have been 
well studied in the literature. 

The layout of this paper is as follows: In section[2]we dis- 
cuss salient details of the simulation including the PNKll 
model and how we have augmented it to allow for black 
hole tracking. In section |3] we compare the behaviour of the 
merger simulations and the sensitivity of final state diag- 
nostics, specifically the black hole mass, as a function of the 
model parameters. We end with a brief review. 



2 NUMERICAL SIMULATIONS 

We test the ADP method in a simulation of a major merger 
of two Milky Way-sized galaxies. Each galaxy contains a 
dark matter halo, hot gas halo, stellar disc, gas disc, a stellar 
bulge, and a black hole. The galaxies are initially separated 
by 70 kpc, and placed on a parabolic trajectory around one 
another similar to the merger considered in [Springel et al.| 
(2005b I. The construction of the galaxies and definitions of 



the structural parameters are given in Appendix [Xj while 
the parameters of each component are presented in Table 
Al| The masses and particle numbers used are summarized 
in Table [H 

The simulations were run using the parallel version of 
Hydra ( Couchman et al.|1995 Thacker & Couchman 2006 1 , 
which uses an Adaptive Particle-Particle, Particle-Mesh al- 
gorithm ( Couchman||1991 1 to calculate gravitational forces 
and the Smooth Particle Hydrodynamics method (SPH; 
Gingold & Monaghan 1977 Lucy 19771 to calculate gas 



forces. The star formation algorithm is described injThacter] 
& Couchman|((2000l. 



To implement an AGN feedback algorithm motivated 
by the PNKll model, we have had to make some modifi- 
cations and additions. We have modified how feedback is 
returned to the local gas, and added both a black hole ad- 
vection algorithm and a black hole merger algorithm. None 
of the issues related to black hole advection or mergers are 
addressed within the PNKll paper since they consider ac- 
cretion of gas on to a black hole that is embedded within an 
initially spherical gas cloud. 
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Fiducial resolution 


Low resolution 






mass 


Particle mass 


Number of 


Particle mass 


Number of 






Mq) 


(105 M0) 


particles 


(105 M0) 


particles 


Dark matter halo 




89.92 


11.75 


765 000 


89.92 


100 000 


Hot gas halo 




0.60 


0.36 


165 343 


2.77 


21 619 


Stellar bulge 




1.34 


2.37 


56 649 


18.10 


7 407 


Stellar disc 




3.56 


2.37 


150 375 


18.10 


19 662 


Gas disc 




0.54 


0.36 


150 375 


2.77 


19 662 


Black hole 




10-^ 


1.00 


1 


1.00 


1 


Total 




95.96 




1 287 743 




168 351 



Table 1. Component breakdown for each galaxy. The accretion disc initially has zero mass, thus the initial mass of the ADP is the same 
as the black hole. 



2.1 Implementing the PNKll model: Approach 
to accretion 

The black hole accretion algorithm is the same as in PNKll. 
In this method the accretion disc particle (ADP) is a coUi- 
sionless sink particle containing both the black hole and its 
tightly bound accretion disc. The mass of the ADP is given 
by 



MadP =AfBH + A/disc, 



(3) 



where AIbh is the mass of the black hole and Afdisc is 
the mass of the accretion disc. As in PNKll, we initialise 
Afdiac = 0, while we choose the initial black hole mass in 
each galaxy to be 10^ Mq. The ADP has an associated 
'smoothing length', /iadp, which is calculated at every it- 
eration by /lADP = max (/lADP, /imin), where a sphere with 
radius 2/iadp around the ADP includes 60 gas particles, and 
hmin is the smallest resolved smoothing length in the SPH 
solver. The smoothing length is used to calculate the gas 
properties at the location of the ADP, which is used for anal- 
ysis and returning feedback energy. The ADP is also given 
an accretion radius, R^cc, which is a free parameter of the 
model that determines at what radii particles are considered 
to have accreted on to the accretion disc. Nominally, and in- 
deed in PNKll, this radius is expected to be on the order of 
a pc. While merger simulations can reach these small radii 
by considering arbitrarily small values of 7?acc, it is worth 
emphasizing that the properties of the gas on these scales is 
not resolved. Indeed this is a generic problem with feedback 
models in general, namely that they use input parameters 
at the edge of model resolution. However, since this is an 
unavoidable problem, we continue on with this issue noted 
and examine the impact of varying the value of J?acc by up 
to a factor of 10. 

Whenever a gas particle comes within Ra.cc of the ADP, 
it is instantly accreted on to the accretion disc, thus Afdisc 
Mdiac + Ttgas, where rrigas is the mass of the gas particle. 
The accreted particle's mass and momentum are added to 
the ADP, and then the accreted particle is removed from 
the simulation. This capture rate is not limited in any way, 
and solely depends on the particles' relative locations. Once 
accreted on to the accretion disc, it takes time for the gas 
to be transported through the disc so that it can finally be 
accreted on to the black hole. This time delay is on order of 
the disc viscous time, tvisc, which must be larger than the 
dynamical time at the accretion radius; PNKll set tvisc as 



a free parameter, but argue that it should be tvisc ~10-100 
Myr in galaxy formation simulations. 

The mass accretion rate on to the black hole from the 
accretion disc is given by, 



A/bh = min ( M^lEl ^ M^^^ 



(4) 



To preserve mass continuity, the mass added to the black 
hole is removed from the accretion disc although the overall 
ADP mass remains the same. As in most numerical imple- 
mentations of AGN feedback, the accretion rate on to the 
black hole is Eddington-limited; this moderates the growth 
of the black hole based upon a physical limit rather than 
just from the numerical accretion on to the disc. 

2.2 Implementing the PNKll model: Changes 
and additions for merger simulations 

2.2.1 Feedback 



PNKll use the same feedback scheme as described in |Nayak^ 
shin & Power (20101. In this scheme, feedback energy is re- 



turned by adding wind particles which are radially directed 
away from the black hole; the wind particles have momen- 
tum Pwind ~ O.lmgaacr, whcre a is the velocity dispersion 
of the host galaxy. In our simulation, it is impractical to 
continually add wind particles since the particle load in the 
solver will climb rapidly, so we instead follow the wind pre- 
scription used in Debuhr et al. (20111. Here, the feedback 
rate is given by 



L 



P 



(5) 



where r = 10 is the assumed infrared optical depth, and 
L = erA^BHC^ is the luminosity, where tr = 0.1 is the ra- 
diative efficiency. The momentum is returned radially, such 
that every gas particle within the ADP's radius of influence, 
Tinf = 2/!,ADP, receives an equal acceleration. We reiterate 
that rinf has no explicit dependence on i?acc and is being 
recalculated at every iteration. 

2.2.2 Black hole advection and mergers 

The black hole advection algorithm is a modified version of 



that presented in Okamoto et al. ( 2008 1 and used in Model 
WT in [Wurster fc Thacker| ( |2013[ ). Here, the ADP is dis- 
placed towards the centre of mass of the sphere with radius 
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rinf which is centred on the AD P. The distance it is displaced 
is 



Al = min(0.10/iADP, 0.30 |v| dt, dcm) 



(6) 



where v is the velocity of the black hole, dt is the time-step, 
and dcM is the distance from the black hole to the centre 
of mass. The coefficients are based upon those in |Okamoto| 
et al. (20081, but modified for our higher resolution. Follow- 
ing the approach of Okamoto et al. (20081, we have selected 



the parameters to suppress the oscillations of the black hole 
particles such that the maximum amplitude of the black 
hole oscillation between the start of the simulation and first 
periapsis is less than 0.25 per cent of the core radius. This 
method continually displaces the black hole towards the cen- 
tre of mass to counter any two-body interactions that may 
try to displace the black hole from the bottom of the poten- 
tial well; limiting the distance preserves the possibility of the 
black hole remaining in a gas void. One option we considered 
but did not implement was to couple the ADP to the gas 
particle that has the lowest potential energy, is within rinf of 
the ADP and satisfies «rci < 0.25cs, where Vrai is the relative 
velocity between the ADP and the gas particle and Cs is the 
local sound speed; generally, this can lead to large artificial 
displacements and prohibits the black hole from remaining 
in a gas void, but with this specific accretion algorithm, it 
would artificially increase the accretion rate by arbitrarily 
moving the ADP within -Race of a gas particle. A second op- 
tion we considered was to use a tracer particle whose mass is 
~100-1000 times greater than any other particle; although 
this produces a smooth black hole trajectory, the additional 
mass affects the evolution time when compared to non-tracer 
particle models. 

Finally, when the ADP's pass within each others 
smoothing lengths and have a relative velocity less than the 
local sound speed, they are assumed to instantly merge. The 
merged black hole has the combined mass of the progenitor 
black holes, and the merged accretion disc has the combined 
mass of the progenitor discs, thus all masses are conserved. 
This is similar to the numerical merger procedure used in 
Springel et al.l (|2005bl). 



2.3 Parameter space and resolution dependence 

To understand the impact of the two free parameters, i?acc 
and tvisc , we ran a suite of simulations within this parameter 
space and considered two separate resolutions for a total 
of 17 simulations. We plot where the simulations lie in the 
parameter space in Fig. [l] where each point corresponds to 
a model. Our parameter space is 

(i) resolution £ {low, fiducial}, 

(ii) i?acc//imin G {0.02,0.05,0.10,0.20}, and 

(iii) tvisc/Myre {1,5,10}. 

In Table [2) we convert the accretion radius to 
physical units for both resolutions. Our model 
naming convention is understood as follows: 
'PNK^e solution 

r(100-Racc/ftmin)t(ivisc/Myr)'. If we do 
not include a resolution when referring to a model, then 
we are referring to both resolutions. In terms of wall-clock 
time, despite having a modest particle content, the fiducial 
resolution models still take over a month to run due to the 
large number of time steps required, the lack of multiple 



^acc 

(Viii) 


Fiducial resolution 
(pc) 


Low resolution 
(pc) 


0.02 


0.731 


1.827 


0.05 


1.827 


4.567 


0.10 


3.654 


9.133 


0.20 


Not run 


18.266 



Table 2. The accretion radii, -Race, in our parameter space, given 
in simulation and physical units. 



Physical • 
Unphysical >K 



Fiducial 



Low 




0.0: 
^acc (^min) 



(Myr) 



Figure 1. The parameter space that we tested. Each point cor- 
responds to a model; the physical/unphysical designation will be 
discussed in ijsj 



time-steps in the Hydra solver and also an 0{n^) slowdown 
that occurs as large numbers of SPH particles reach the 
minimum smoothing length of the solver. We thus have a 
limited number of fiducial resolution simulations. 



3 RESULTS 

3.1 General evolution 

Each of our models was evolved for 1.5 Gyr through a merger 
event, similar to that of [Springel et ah] ( |2005b[ ). We found 
that each model followed a similar qualitative history; in 
Fig. [2] we plot the evolution of the gas column density of 
Model PNKfr05t05. From this evolution, we see four signif- 
icant epochs; first periapsis, apoapsis, second periapsis and 
core merger, occurring at ~170, 480, 880 and 990 Myr, re- 
spectively; the low resolution models reach second periapsis 
and core merger 25 and 6 Myr, respectively, earlier than their 
fiducial resolution counterparts. In all our fiducial resolution 
models, the black holes merge at 1.01 Gyr, and their low 
resolution counterparts merge 10-60 Myr later. The later 
merger in the low resolution models is a result of the binary 
black holes oscillating about one another with a greater am- 
plitude than in the fiducial resolution models; this is from 
both the low resolution black holes being more massive by 
this epoch (see section [3^ , and there being fewer particles 
to induce drag on the black hole system. 
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Figure 2. Evolution of gas column density for Model PNKfr05t05. Times from the onset of the simulation are listed in each frame; each 
frame measures 100 kpc per side, with an image resolution of 98 pc/pixel. 



The final remnant is a reformed gas disc and flattened 
stellar ellipsoid surrounded by a hot gas halo; the remnant 
is further discussed in section |3.7| The final gas discs in 
the fiducial resolution models have scale lengths of approx- 
imately 0.5 kpc, compared to the initial scale length of 2.46 
kpc and are essentially an order of magnitude less massive 
than the discs in the initial conditions. Of the gas that was 
initially in a disc, 89 per cent of it is either accreted on to the 
black hole or converted into stars. With a final star forma- 
tion rate of less than 0.5 MQ/yr, the final configuration is a 
red and dead elliptical, as expected ( [Springel et al.|2005a[ ). 

To quantify this evolution, we have plotted the total 
black hole mass, the accretion rates on to the black hole, and 
the gas density and gas temperature within rinf of the ADP 
for Model PNKr05r05 in Fig. |3] the accretion rates and gas 
properties are geometrically averaged over both black holes 
and plotted in bins of 10 Myr. 



For the first ~300 Myr, rinf > 2/ii„in, thus there are 
only 60 gas particles within this radius; after ~300 Myr, 
rinf = 2/iniin, and there are several thousand particles within 
this radius. There is an increasing accretion rate after first 
periapsis, which peaks near apoapsis; there is a second peak 
at core merger. The gas density in the core also increases 
from first periapsis to apoapsis, and again at core merger. 
The core temperature stays relatively moderated after first 
periapsis. The initial decline in gas density and temperature 
is a result of the feedback energy carving out a void around 
the black hole in the initial density field which is quite sym- 
metric. The SPH particles that contribute to the density 
and temperature at the location of the black hole are at 
r ~ 0.9rinf, thus provide only a minimal contribution. Thus, 
these low densities and temperatures are essentially a result 
of the smoothing kernel breaking down due to a poor dis- 
tribution of the particles within 2/iadp. At later times the 
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Fiducial Resolution 



Low Resolution 




a VISCOUS accretion rate, 



0.2 0.4 0.6 0.8 1 1.2 1.4 0.2 04 0.6 0.8 1 1.2 1.4 
Time (Gyr) Time (Gyr) 

Figure 3. Top to bottom: Total black hole mass, accretion rates 
on to the black hole, gas density and gas temperature within 
rjnf of the ADP/black hole particle. The bottom three rows are 
geometrically averaged over both black holes in bins of 10 Myr. 
The left panels are from the fiducial resolution simulations and 
the right panels are from their low resolution counterparts. The 
vertical line of the same lincstyle indicates when the black holes 
merge. 



number of particles increases and symmetry is lost which 
prevents any spurious temperature and density values from 
being calculated. 



3.2 Comparison to other models in the literature 

We briefly compare Model PNKr05t05 to two other models 
found in the literature: Model SDH, which is designed to 
reproduce the algorithm used in [Springel et al.| ( |2005b[ ) and 
Model DQM, which is designed to reproduce the algorithm 
used in iDebuhr et al.l (|2011|). Model SDH uses the Bondi 



(20111. 

accretion rate given in (1|, where a = 100, to calculate the 
rate of gas accretion on to the black hole; half a per cent of 
the accreted mass is converted in to feedback energy which 
is returned thermally by adding energy to the nearby gas 
particles, following the smoothing kernel. Model DQM uses 



Afvisc = 3tt5T, 



(7) 



where 5 — 0.05 is the dimensionless viscosity, E is the mean 
gas surface density, and O. = ^ GM/r'f^f is the rotational 
angular velocity of the gas, which is based upon the multi- 



scale SPH simulations by Hopkins & Quataert (20101; the 



feedback algorithm is same as in our PNK models. We plot 
the total black hole mass, the accretion rates on to the black 
hole, and the gas density and gas temperature within rinf of 
the black hole for Models SDH and DQM in Fig. [S] 

We briefly compare Model PNKr05t05 to two other 
models found in the literature: Model SDH, which is de- 
signed to reproduce the algorithm used in [Springel et al.| 
( |2005b| ) and Model DQM, which is designed to reproduce 
the algorithm used in Debuhr et al. (2011 1. Model SDH uses 



the Bondi accretion rate given in ([l| to calculate the rate 
of gas accretion on to the black hole; the feedback energy 
is returned thermally by adding energy to the nearby gas 
particles. Model DQM uses a viscous accretion rate, which 
is based upon the multi-scale SPH simulations by |Hopkins| 
& Quataert ( 2010 1; the feedback algorithm is same as in our 



PNK models. We plot the total black hole mass, the accre- 
tion rates on to the black hole, and the gas density and gas 
temperature within rini of the black hole for Models SDH 
and DQM in Fig. [S] 

Although all three accretion and feedback algorithms 
make different physical assumptions, the final total black 
hole masses in the fiducial (low) resolutions are equal within 
a factor of 1.52 (4.62). The accretion histories, however, vary 
considerably. As expected. Models SDH and DQM have con- 
tinual accretion which starts immediately, but the accretion 
rate in PNKr05t05 is punctuated by periods of nearly zero 
accretion, especially at early times. The latter is a result of 
feedback from the first few accretion events creating a small 
and transient void around the ADP, as can be seen in the 
core density profile. After first periapsis, tidal torques are 
large enough and feedback energy is low enough to allow 
additional accretion events. 

The core density of Model DQM is the lowest of these 
three models; in Model DQMf, kinetic feedback builds up 
in the gas, and after 60 Myr, a void of ~ 0.55rinf is formed 
and maintained for the duration of the simulation. Although 
Model PNKr05t05 uses the same feedback algorithm, the 
initially lower accretion rate prevents this build up of ki- 
netic feedback, thus this void never forms, hence the higher 
core density. The core temperature is highest for Model 
SDH, where thermal feedback energy directly increases the 
temperature of the gas. The core heating of Models PNK 
and DQM is from shock heating and star formation; the 
higher accretion and star formation (not shown) rates in 
PNKr05t05 after first periapsis results in more heating than 
in Model DQM. 



3.3 Acceptable parameter ranges 

We ran 12 low resolution simulations to probe our entire 
J?acc-ivisc parameter space. At our fiducial resolution we 
only tested five models due the time required to complete 
individual simulations. After analysing the results, each sim- 
ulation was classified as either physical or unphysical, as is 
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indicated by the dot type in Fig. [T] The basis for our def- 
inition of physical or unphysical rehes upon a combination 
of structural evolution and how close the final remnant lies 
to the Mbh-o" relationship. We note that relying upon the 
Mbh-o" relationship is perhaps not a strong constraint be- 
cause the masses of the seed black holes could be changed, 
however it remains a commonly used approach in these simu- 
lations and we thus proceed cautiously with this issue noted. 

For -Race values that are very small the accretion rate on 
to the black hole will be limited. The expectation in this case 
is that the resulting remnant will fall below the A/bh-o" rela- 
tion. Our numerical experiments show this to be true for all 
the models with -Race ~ 0.02/iniin, with the black hole mass 
in PNKfr02t05 being 4.5 times lower than expected, assum- 
ing the velocity dispersion remains constant. These mod- 
els also suffer from sensitivity to operation ordering in the 
solver. At a given resolution, we would expect the absorp- 
tion of particles to occur at the same time provided that the 
operations in the solver are executed in the same order, and 
hence produce the same round-off. However, when running 
in parallel with dynamic load-balancing this is no longer the 
case and round-off in accumulations can produce subtle dif- 
ferences around or just above machine precision. Thus, with 
too small of an accretion radius, pseudo-random issues due 
to machine precision dominate the accretion, and the results 
of a given simulation are not easily reproducible. On the ba- 
sis of lying beneath the Mbh-ct relation and the exhibited 
numerical sensitivity, we classify all of our -Race = 0.02/imin 
models as 'unphysical'. 

At the other extreme of our parameter space, there can 
be unphysical results if the accretion rate on to the disc 
is too high. If -Race is too large, then many particles can 
pass within this radius and be accreted on to the disc, re- 
gardless of the amount of feedback. Likewise, a large t^isc 
means that the feedback energy is being returned at a lower 
rate, thus particles have less of an obstacle to overcome to 
pass within 7?aec when influenced by outside forces (i.e. tidal 
torques from the interacting galaxies) . Both result in a large 
accretion rate on to the disc, where the accretion disc mass 
can surpass the mass of the black hole (see bottom right 
panel of Fig. [4] the remaining three panels show the profiles 
for physical models, where the peak accretion disc mass is 
~ 5 X IO^Mq < Mbh). This unreasonably large accretion 
rate on to the disc yields a large and continual accretion rate 
from the disc on to the black hole, which results in a vast 
amount of feedback energy. Once a critical amount of kinetic 
feedback is injected into the gas, namely that sufficient to 
overcome the nuclear gravitational binding energy, a large, 
unphysical void is formed in the galaxy. This in turn sup- 
presses the accretion rate on to the ADP, but not from the 
disc on to the black hole. The accretion disc is ultimately 
depleted, the feedback ends and the void recoUapses to begin 
another cycle. In situations where a large radius and large 
viscous time-scale are included it is possible to remove the 
gas from the system. 

These unphysical voids often foreshadow the total dis- 
ruption of the remnant. Starting at 1.04 Gyr in Model 
PNKfrOStlO, there is a catastrophic explosion, and the sys- 
tem is totally disrupted within 60 Myr. In other models, 
specifically our low resolution models, the core merger causes 
an increase in the accretion rate, leading to a final feedback 
episode that disrupts the system, although often less catas- 



PNK^05110 (unphysical) 
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Figure 4. Total black hole mass (solid red line) and total ac- 
cretion disc mass (dashed blue line; geometrically averaged over 
both black holes and plotted in bins of 10 Myr) for four fiducial 
resolution models. The vertical line indicates when the black holes 
merge. 



trophically than Model PNKfr05tlO. While the system is 
highly disrupted, there are often no noticeable voids, and 
the process is 'gentle' enough that the gas begins to recon- 
dense into a disc and cloud. We cautiously and liberally 
classify these results as physical. 

In Fig. [5] we have plotted the actual and the Edding- 
ton accretion rates for six models; we also give the length of 
time each system is undergoing Eddington-limited accretion 
in the top left corner of each panel. The unphysical model 
PNKfrOStlO is undergoing Eddington-limited accretion for 
11.7 per cent of the time, thus a (relatively) high percentage 
is a possible indicator of an unphysical model. However, the 
unphysical Model PNKfr02t05 undergoes Eddington-limited 
accretion 6.3 per cent of the time, which is higher than the 
per cent of time the physical Model DQMf is undergoing 
Eddington-limited accretion. While there is a correlation be- 
tween the amount of time spent in Eddington-limited accre- 
tion and the final black hole mass for the PNK models, this 
is not true when comparing across different models. Further, 
these numbers indicate that relying upon the Mbh-ct rela- 
tionship appears to place the amount of time in Eddington- 
limited accretion in a comparatively narrow band. It is thus 
difficult for us to draw detailed conclusions about using the 
amount of time spent in Eddington-limited accretion as a 
strong discriminant between the physical or unphysical na- 
ture of models. 

Ultimately, the amount of time accreting at the Ed- 
dington limit is very closely related to the duty cycle of the 
quasar. In the context of models such as [Small fc Bland-] 
|ford| ( |1992[ ), where the population of SMBH is predicted as 
a function of luminosity (or mass) and time, the duty cycle 
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T(Mj,jj) = 0.094 Gyr 



PNK,r05tlO (unphysical) 
T(ME<u) = ai75 Gyr 




Fiducial Resolution 



Low Resolution 



DQM, 

T(MEjd) = 0.059 Gyr 





SDH, 

T(MEjd) = 0.103 Gyr 
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Time (Gyr) Time (Gyr) 

Figure 5. Actual (solid red line) and Eddington (dashed blue 
line) accretion rates for six models, geometrically averaged over 
both black holes and plotted in bins of 10 Myr. The vertical line 
indicates when the black holes merge. Model names and length 
of time the system is undergoing Eddington accretion are listed 
in the top left of each panel. 




0.2 0.4 0.6 0.8 1 1.2 1.4 02 0.4 0.6 0.8 1 1.2 1.4 
Time (Gyr) Time (Gyr) 

Figure 6. The total black hole mass, accretion rates on to the 
black hole, gas density and gas temperature within r^^f of the 
ADP, as in Fig. [3] we fix tvisc = 5 Myr and vary Race- 



directly impacts the global population because the luminos- 
ity function is a product of the number density, mass and 
duty cycle. Models of this type have been extensively devel- 



oped by Shankar (e.g. Shankar et al. 2013 Shankar et al 
|2010[ ) and show that while Eddington accretion is likely to 
be more common at lower masses (Af < 10^'^ M^0)i once the 
SMBH reaches 10* Mq, matching the luminosity function 
requires that growth at the Eddington rate must be con- 
siderably less frequent or radiative efficiencies must be low. 
While obviously not directly equivalent, matching the A4"bh- 
cr relationship appears to place a very similar constraint on 
the allowed duty cycle. 

Lastly, the gas discs in the simulation are quite stable 
and in the absence of any additional supply of infalling gas 
actually decrease their mass during the simulation by 63 
per cent prior to the start of the main merger. The stabil- 
ity of the discs is a product both of our choosing a com- 
paratively stable initial configuration and also that the SN 
feedback routine keeps the gas comparatively well supported 
against cold collapse. Thus the discs do not fuel the black 



hole through large scale fragmentation due to them becom- 
ing unstable. 



3.4 Parameter sensitivity: i?acc 

We have four models (two models each of fiducial and low 
resolution) with tvisc ~ 5 Myr that produce physically plau- 
sible results. In Fig. [6] we have plotted total black hole mass, 
the accretion rates on to the black hole, and the gas density 
and gas temperature within rinf of the ADP. As expected, 
the first accretion event occurs in Model PNKrl0t05, and 
Model PNKrlOtOS has a slightly higher accretion rate for 
most of the simulation. By increasing the radius of accretion 
by a factor of two (hence the volume in which a particle can 
be accreted, T4cc, by a factor of eight), there is only a fac- 
tor of 1.35 and 1.70 increase in final total black hole mass 
for fiducial and low resolutions, respectively. The gas den- 
sity and temperature behaviour within rinf are very similar 
at both resolutions. The value of rinf has no explicit depen- 
dence on i?acc and all models produce the same qualitative 
behaviour for rinf, which settles to ~ 2/iinin between first 
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Figure 8. The gas accretion rate on to one black liole for ttiree 



low resolution models, holding _Ra 



0.05/ln 



We show the 



initial accretion event (top) and the accretion just prior to core 
merger (bottom). 



Figure 7. The total black hole mass, accretion rates on to the 
black hole, gas density and gas temperature within rj^f of the 
ADP, as in Fig. [3] we fix -Race = O.OS/imin and vary tvisc- 

periapsis and apoapsis. We find that -Race has neghgible in- 
fluence on rinf and a minimal impact on the final black hole 
masses. Thus, we conclude that, as long as the parameters 
are in the physically acceptable range, results are not overly 
sensitive to the exact value of -Race- 

3.5 Parameter sensitivity: tvisc 

We have five models with -Race ~ 0.05/imin (two at the fidu- 
cial resolution and three at low resolution). In Fig. [7j we 
have again plotted total black hole mass, the accretion rates 
on to the black hole, and the gas density and gas temper- 
ature within Tinf of the ADP. The total black hole masses 
in the fiducial resolution simulations are similar through- 
out the simulations, and their final masses differ by 13 per 
cent. At low resolution, the total black hole masses for mod- 
els PNKir05t05 and PNKirOStOl match within 18 per cent 
at the end of the simulation; if we include PNKirOStlO (a 
model that produces unphysical behaviour at the fiducial 
resolution), then the low resolution black hole masses differ 
at most by a factor of 1.68. 



In Fig. [S] we plot two segments of the accretion rate 
for one black hole in our three low resolution models with 
Ra.cc ~ 0.05/imin. As expected, the first accretion event hap- 
pens simultaneously in each model (top panel). This is fol- 
lowed by a period of Eddington accretion, which stops first 
for tvisc ~ 10 Myr and last for tvisc = 1 Myr depending 
upon when the -Mdisc/ivisc rate falls below the Eddington 
rate. Once the Eddington accretion has stopped, the rate of 
feedback (where p oc M) differs for each model, but at least 
at this early time the accretion rate on to the black hole 
can still be calculated: an exponentially decaying rate with 
the time dependence given by t/tvisc- Although PNKr05t01 
accretes at the Eddington limit very slightly longer than the 
other two models, the faster decrease in the accretion rate 
means that feedback effectively stops sooner, thus giving 
particles more time to slow down and fall within -Race- 

At late times, each accretion event does not add enough 
mass to the disc to permit Eddington accretion (bottom 
panel); however, there is a local accretion peak followed by 
exponential decay. As previously discussed, the faster de- 
cay of tvisc = 1 Myr allows for rapid accretion events on to 
the disc; this creates large variances in the accretion rates, 
spanning a few orders of magnitude. The slower decay of 
tvisc ~ 10 Myr yields a moderated accretion rate, thus there 
is only a factor of a few between the local peak accretion rate 
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and the minimum rate prior to the next accretion event. The 
time-average accretion rate yields a higher accretion rate for 
the lower viscous time-scales, resulting in more feedback be- 
ing returned to the gas which can, at least temporarily, expel 
considerable amounts of gas from the system. Summarily, 
we find a general trend to more massive black holes with 
decreasing tvisc- 

The torques from the interacting galaxies modify this 
argument by introducing disc instabilities Although the 
feedback rate of the tvisc = 10 Myr models can hinder gas 
from accreting on to the disc during a quiescent phase, there 
is not enough cumulative feedback to prevent a large gas flow 
on to the disc produced by the disc instabilities. This results 
in a large accretion epoch, followed by powerful outbursts 
of feedback energy. Thus, the tidal torques at first periapsis 
and core merger are strong enough to over come the feedback 
and cause a catastrophic (or near catastrophic) outburst, as 
in Model PNKfrOStlO. 

Similar to our discussion in section [374} we find quanti- 
tatively similar gas behaviour within rinf, which is a result 
of the rinf's being similar across all models. As with i?acc, 
we conclude that, as long as fvisc lies with the physically 
acceptable range, then the simulations are comparatively 
insensitive to its exact value (at least to within 1/2 an order 
of magnitude at the fiducial resolution). 

3.6 Resolution sensitivity 

We have three models that are deemed to have physically 
plausible results at both resolutions. In Fig. [9j we have plot- 
ted total black hole mass, the accretion rates on to the black 
hole, and the gas density and gas temperature within Vn-ii of 
the ADP. When we directly compare two models at different 
resolutions, there are significantly more differences in accre- 
tion and evolution histories than we do for when the model 
parameters are changed (at least for the acceptable range). 
As would be expected, in all cases the fiducial resolution 
model accretes the first particle, while the lower resolution 
simulations produce a larger jump in the accretion rate at 
early times due to the particles in those simulations being 
more massive. By the time the galaxies reach first periapsis 
at 166 Myr, the low resolution gas has a greater kinetic feed- 
back obstacle to overcome before it can pass within -Race of 
the ADP. Thus, the first low resolution accretion event im- 
mediately begins to moderate the accretion on to the disc, 
and this moderation continues to persist after first periapsis. 

Because it has a smaller kinetic feedback obstacle to 
overcome, the fiducial resolution gas can more easily fall 
within -Race of the ADP after first periapsis, leading to a 
large and essentially unmoderated accretion on to the disc. 
As the gas accretes from the disc on to the black hole, it mod- 
ifies the environment to suppress further accretion events 
on to the disc; however, the disc remains massive from the 
previous accretion episode. Thus, there is a decrease in gas 
density in the core, but the accretion on to the black hole 
continues as the gas in the disc is continually transferred to 
the black hole. This major epoch of accretion on to the disc 
after first periapsis results in a steeper black hole growth be- 
tween first periapsis and apoapsis for the fiducial resolution 
models. 

For all three models, the fiducial resolution models have 
the expected higher core gas densities. Since the feedback is 
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Figure 10. Our numerical Mbu^ct relation, along with the ob- 
served relation (red solid line) and the one-sigma scatter (green 
dashed line) from |Giiltekin et al.|(2009| l. For six selected models, 
the solid dot/triangle represents the average cr of 1000 random 
lines of sight, the horizontal bars represent the range of all calcu- 
lated cr's, and the remaining three symbols on the horizontal line 
represent cr taken along the itx— , iy— and ±2— lines of sight. For 
the remainder of our physical PNK models, we have only plotted 
the average cr of 1000 random lines of sight. We have explicitly 
labelled which points/set of points corresponds to which models. 



being returned kinetically, shock heating and star formation 
are the primary heating mechanisms, and both scale with 
resolution in turn producing similar core temperatures. 

Thus, based upon the two resolutions we test, we can 
conclude that resolution has a greater impact on the results 
than the values of the free parameters, -Race and tvisc. How- 
ever, these differences are predictable since there are more 
epochs of discrete behaviour in the fiducial resolution mod- 
els than in the low resolution models, allowing for a more 
continuous modelling of the evolution. Also, each fiducial- 
low resolution pair is more similar to one another than to 
Models SDH or DQM of the same resolution, indicating that 
the model can be distinguished from other models even at 
low resolution. 



3.7 Final states 

3. 7. 1 Stellar remnant 

A common test of numerical accretion and feedback algo- 
rithms is the -Mbh-ct relation. For all of our physical models, 
we have calculated the stellar velocity dispersion around ev- 
ery black hole using the method described in |Debuhr et al.| 
(20111. These velocity dispersions are averaged over 1000 
random lines of sight, and are plotted on the -Mbh-o" rela- 
tion in Fig. |10[ For the six models in Fig. [9] we also plot 



the full range of velocity dispersions, as well as those taken 
preferentially along the ±x-, ±y- and ±2:-lines of sight. The 
large range of cr is a result of the highly triaxial stellar rem- 
nant; see Fig. |11[ where we have plotted a fiducial and low 
resolution stellar remnant. The fiducial resolution models 
have average values very near the observed Mbh-o" relation, 
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Figure 9. The total black hole mass, accretion rates on to the black hole, gas density and gas temperature within rj^f of the ADP, as 
in Fig.|3] Each column plots data for a different set of iJacc and ivisc- 
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Figure 11. Stellar column density of the remnants of PNKfr05t05 (left) and PNKir05t05 (right), taken at 1.5 Gyr. All stellar remnants 
(including Models DQM and SDH) of each resolution are similar, with the low resolution remnants yielding more elliptical bulges than 
their fiducial resolution counterparts. Top Row: Face-on view with each frame measuring 100 kpc per side, with image resolution of 195 
pc/pixel (left) and 390 pc/pixel (right). Bottom: Face-on and edge-on view of both models, with each frame measuring 20 kpc per side, 
and image resolution of 39 pc/pixel (78 pc/pixel) for the fiducial (low) resolution models. 



and a range that nearly stays within the one-sigma scat- 
ter. The low resolution stellar remnants are more elliptical 
than their fiducial resolution counterparts, thus they have 
a greater range of a's. While some of the average cr's fall 
outside of the one-sigma range, they all obtain a's near the 
observed Mbh-o" relation if the line of sight is near the plane 
of the disc. To verify these results, we have recalculated the 
velocity dispersions of the fiducial resolution models based 
upon the gravitational potential of the black hole, and found 
these values to be consistent within five per cent of the av- 
erage values reported above. 



Model 



Disc mass (10* Mq) 



PNKfr05t05 

PNKfr05t01 

PNKfrl0t05 

DQMf 

SDHf 



4.18 
3.93 
3.31 
6.42 
2.70 



Table 3. The gas mass of the reformed discs at 1.5 Gyr for five 
fiducial resolution models. 



3.7.2 Gas properties 

In Figs. [12] and |13[ we plot the gas column density of ten 
remnants at 1.5 Gyr. These show, respectively, the inner 100 
and 20 kpc of six PNK models (three models as each resolu- 
tion) and Models DQM and SDH. As discussed in sections 
|3.4| and |3.5| there is much qualitative similarity amongst the 
three fiducial remnants of the PNK models, as well as Mod- 
els DQMf and SDHf. In all cases, there is a condensing gas 
cloud with a reformed disc. 

The radius, surface density profiles and total gas mass 
of the reformed discs in the PNK models are only slightly 
dependent on the free parameters; the surface density profile 
is plotted in the top panel of Fig. [14] and the gas disc masses 
are given in Table [3] 



After core merger, PNKfrlOtOS has a slightly larger ac- 
cretion rate than PNKfr05t05 due to its larger accretion 
radius, this leads to more feedback energy and a slightly ex- 
tended disc with lower surface density. Model PNKfrOStOl 
has a higher accretion rate at core merger than PNKfr05t05, 
which results in a small outburst at 1 Gyr, expelling some 



gas from the system. By 1.5 Gyr, the remaining gas is less 
bound leading to the slightly larger radius and lower surface 
density. 

The fiducial resolution PNK models have different sur- 
face density profiles than Models DQMf and SDHf; see the 
bottom panel of Fig. [l4] and note the difference range on the 
vertical axis. The disc in Model DQMf has a dense core and 
a moderately dense torus (resulting in a high-mass disc), 
whereas the profile for Model SDHf is a torus due to low 
level feedback activity carving out a region of the disc. 

The three fiducial resolution PNK remnants are quali- 
tatively more similar to one another than to their low reso- 
lution counterparts. Both PNKir05t01 and PNKir05t05 un- 
derwent a major outburst starting at 1.05 Gyr, blowing away 
much of the gas and preventing the reformation of the gas 
disc. In the subsequent few 100 Myrs, the gas begins to re- 
condense into the cloud presented here. In PNKirl0t05, this 
major outburst never occurs, thus allowing the disc to re- 
form. 



ADP Accretion in Major Merger Simulations 13 



PNKfl 
• 


05t05 


PNK|r05t05 
* « 


PNKfl 


r05t01 


PNK|r05t01 


T 








PNKfl 


rl0t05 


PNK|rl0t05 

♦ 


DQMf 

4 

• 


DQ 


M, 


SDHf 

b — 


SDH| 

• 


PL 



h-» o O M ho w 

(n cji ui ui cn CJ1 



log ( Column Density qas/(M sun 'pc^) ) 



Figure 13. Gas column density of ten remnant, taken at 1.5 Gyr. Each pair contains a face-on and edge-on view of the central 20 kpc. 
The image resolutions are 39 pc/pixel (78 pc/pixel) for the fiducial (low) resolution models. 
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Figure 12. Gas column density often remnant, taken at 1.5 Gyr. 
Eacti frame measures 100 kpc per side, with image resolution 
of 195 pc/pixel (390 pc/pixel) for tiie fiducial (low) resolution 
models. 
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Figure 14. Gas surface density profile, averaged over all az- 
imuthal angles, for the reformed disc of three fiducial resolution 
remnants. The profile is truncated at the edge of the disc, and 
plotted in bins of 49 pc. 



4 CONCLUSION 

We have implemented the accretion disc particle (ADP) 



method of Power et al. (20111 into a major merger simu- 



lation of two Milky Way-sized galaxies. We ran five fiducial 
resolution simulations and twelve low resolution simulations 
varying the free parameters i?acc and fvisc- Our primary con- 
clusions are as follows: 

(i) For accretion radii that are too small (i.e. i?acc ~ 
0.02/iniin), the final black hole mass is far smaller than 
predicted from the Mbh-o" relationship. Thus all of our 
^acc ~ 0.02/iniiii models are classified as unphysical. 

(ii) For accretion radii that are too large (i.e. J?acc ~ 
0.20/imin), the accretion rate on to the accretion disc (hence 
on to the black hole) is unreasonably large. The resulting 
feedback is enough to catastrophically disrupt the system. 
We thus classify all of our iiacc = 0.20/imin models as un- 
physical. 

(iii) For large viscous time-scales (i.e. tvisc = 10 Myr), 
feedback from an accretion event persists long enough to 
hinder secular accretion. Tidal forces from the interacting 
galaxies can overcome the low amount of feedback and fun- 
nel considerable amounts of gas on to the accretion disc. 
Depending on resolution and accretion radius, this short ac- 
cretion epoch can be large enough such that its resulting 
feedback energy can catastrophically disrupt the system. 
Thus a few of our tvisc = 10 Myr models are classified as 
unphysical. 

(iv) The values of -Race and tvisc, assuming they were in 
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the allowed parameter space, had minimal affect on the gas 
properties within rinf of the black hole. This result was ex- 
pected since rinf had no explicit dependence on /?acc or tvisc- 

(v) The exact value of -Race has only a minimal affect 
on the resulting system, assuming that it is in the allowed 
parameter space. By doubling the accretion radius from 
0.05/imin to 0.10/imin, the fiual black hole mass increases by 
only a factor of 1.35 (1.70) for our fiducial (low) resolution 
simulations. 

(vi) As we decrease the value of t^isc the final black hole 
mass increases; the final range of black hole masses spans a 
factor of 1.13 (1.68) for our fiducial (low) resolution models. 
Thus, the exact value of ivisc has only a minimal affect on 
the resulting (physical) system. 

(vii) Decreasing the resolution increases the final black 
hole mass; for any given fiducial-low resolution pair, the final 
black hole mass differs by at most a factor of 1.90. The fidu- 
cial resolution models experience steeper black hole growth 
between first periapsis and apoapsis, and the lower resolu- 
tion models are more prone to major outburst events shortly 
after core merger. Given the parameters tested here, reso- 
lution has the largest impact on the outcome of the model. 
We understand that these differences are unavoidable and 
are a result of a single accretion event in the low resolution 
models returning 7.6 times more energy to the gas than the 
fiducial resolution models. 

In AGN feedback models where the accretion rate is 
dependent on the gas properties near the black hole, the ac- 
cretion and energy feedback rates are never zero, and span 
only a few orders of magnitude; in our simulations, the ac- 
cretion and feedback begins immediately. In Model DQM, 
an initially high accretion rate leads to an immediate and 
dramatic modification of the black hole's environment. In 
Model SDHf, the accretion rate does not become sigiuficant 
until apoapsis, at which point the black hole undergoes a 
major accretion epoch, with its mass increasing by a factor 
of 20.5 in 470 Myr. 

With the ADP model, the initial accretion rate is zero 
or very small. The minimal modification of the ADP's envi- 
ronment by first periapsis allows a large accretion event on 
to the disc, which results in a large accretion rate on to the 
black hole, leading to a very rapid black hole growth; in the 
case of Model PNKfr05t05, the black hole mass increases by 
a factor of 19.6 in 200 Myr. Although these three models, 
Models SDH, DQM and PNK, have very different algorithms 
and different evolution histories, there are similar remnants 
in all three cases. 

The ADP algorithm decouples the accretion rate on to 
the black hole from the gas properties around it and allows a 
high accretion rate even after feedback has prevented addi- 
tional accretion events on to the disc. Yet to have confidence 
in this conceptual model a precise observational strategy for 
determining parameters is necessary (essentially looking at 
the correlation, or lack of, between activity and the nuclear 
gas environment). Recent observational work (Wild et al 



2010P has unveiled how star formation activity and black 
hole growth appear tied together in spheroids. These results 
undoubtedly give useful hints on how black hole growth pro- 
ceeds outside of rapid accretion phases, but fueling of black 
hole via cold gas accretion remains the biggest uncertainty 
in our models currently. To date, studies of nuclear CO mor- 
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Table Al. The chosen parameters for our model galaxies. All 
parameters are defined in Appendix [A] 



phology such as the NUGA project (e.g. Garcfa-Burillo et al. 
|2005[ [van der Laan et al.|201l] ) have found no obvious mor- 
phological links between local AGN activity and mid-scale 
CO morphology. But such studies are obviously limited by 
resolution concerns, as well as probing an entirely different 
AGN luminosity regime. However, ALMA in its full configu- 
ration with the largest baselines will provide resolutions that 
are sufficient to study the nuclear regions the nearest active 
quasar systems. The prospect of examining morphology de- 
pendence in these systems is immensely exciting and should 
provide great insight into their evolution and how they can 
best be modelled. 
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APPENDIX A: GALAXY MODELS 

We construct our model galaxy to have similar mass compo- 



nents to the galaxies presented in Springel et al. ( 2005b I. We 
begin by using the GalactlCs package ( Kuijken & Dubinski 
[Tgisl lWidrow & Dubinski 2005| rWidrow et al.|2008| ) to cre- 
ate a Milky Way-sized galaxy that consists of a stellar bulge, 
stellar disc, and a dark matter halo; this is done through an 
iterative process to produce a self consistent system. The 
free parameters, and the values we chose, are listed in Table 
lATI 



16 J. Wurster and R.J. Thacker 



The dark matter halo profile follows 

C(r;rh,(5rh), (Al) 



Ph = 



1 



4TTal (r/ah)^(l + r/ah)»-T 



where ah is the halo scale length, rh is the cutofi^ radius, 
7 is the central cusp strength and (Th is a (line of sight) 
velocity scale that sets the mass of the halo. The truncation 
function, C(r; rh, 5rh) — |erfc ^ j i smoothly goes from 
one to zero at r = rh over width Srh. 

The stellar disc has a truncated density profile that falls 
off approximately exponentially in R and follows sech^ in z; 
the disc has radial and vertical scale heights -Rd and Zd and 
truncation distances -Rtrunc and ztrunc- The radial velocity 
dispersion is given by cr^(-R) = o-'^Qe"^^^" , where (trq is 
the central velocity dispersion and Re = -Rd for simplicity. 

The stellar bulge density profile is given by 



Pb(r) = pb 



Rc 



(A2) 



which yields a Sersic law for the projected density pro- 
file if p = 1 - 0.6097/n + 0.05563/n^ where n is a 
free parameter. The constant pb is defined using ab = 

|47™6"(''-2'r [n(2 - p)] Rlp^Y^^ , where ol is the depth of 

the gravitational potential associated with the bulge and Rc 
is the radial scale parameter. The variable b is adjusted such 
that _Re encloses half the total projected light or mass. 

The galaxy is then modified in three ways. First, ten 
per cent of the total stellar mass is converted into gas to 
create a gas disc. The gas disc is the same as the stellar disc 
except that it has been reflected in the x = y plane to avoid 
coincidence with the star particles. The gas scale height is 
initially larger than physically motivated, however, cooling 
allows the gas to collapse into a thin disc within a few 10 
Myr. This vertical collapse produces a short transient evo- 
lution of the gas accompanied by a brief increase in the star 
formation rate. At resolutions higher than presented here, 
this produces a strong ring-shaped shock which propagates 
outwards; one possible solution is to reduce the scale height 
of the gas disc, but studies in Williamson & Thacker ( 2012[) 



show that this is not necessary for the resolutions presented 
here. 

Second, a hot gas halo is added following the observa- 



tionally motivated /3-profile (e.g.|Cavaliere k, Fusco-Femiano| 



1976) 



1 + 



(A3) 



where po is the central density, rc is the core radius, and /3 
is the outer slope parameter; we choose rc = 1.75 kpc and 
/3 = 2/3 as done in Moster et al. (2011 1. We set po by choos- 



ing the mass of the hot gas halo within 40 kpc to be equal 



to two per cent of the total disc mass (Rasmussen et al 
2009). To conserve total halo mass, we reduce the mass of 



the dark matter halo by the total mass of the hot gas halo. 
By assuming isotropy and hydrostatic equilibrium, the tem- 



perature profile of the hot gas halo is given by (Kaufmann 
et al.|2007| ) 



T(r) 



/imp 



Phalogas(r) 



phaio gas 



(O^dr, 



(A4) 



where p is the mean molecular weight, rrip is the proton 
mass, /cb is the Boltzmann constants, and M{r) is the total 
mass interior to r. The halo gas is given an initial angular 
momentum which scales with the circular velocity, j(-R) oc 
Rvciic{R), where R is the distance from the spin axis of the 
galaxy ( [Moster et al.|2011 ). 

Finally, a black hole (sink) particle is placed at the cen- 
tre of the galaxy. The galaxy has a total mass of 9.60 x 10^^ 
Mq, 1 287 743 (168 351) particles for the fiducial (low) 
resolution simulations, and a Plummer softening length of 
epi 



ummcr 



= 120 pC (epiummcr = 300 pc) . 
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